function [x,output]=agent12solverho(x0,index,ewsb,waut);
global P beta n m S state agent gammagrid  Y scalefactor lowerbound rho table  ewdev statedev12 statedev13 statedev23
objval = [];
output=[];
options = optimset('LargeScale','off','Display','iter','TolX',0.0001);
[x,fval,exitflag]=fsolve(@objfun,x0,options);
output(6)=exitflag;
function f = objfun(x);
i=index(1);
j=index(2);

consu1 = x(1);
consu2 = x(2);
gammarfun=[x(3) x(4)];

wint=interpolation(scalefactor*ewsb, index,gammarfun);


output(1)=utility(rho,consu1)+beta*wint(1)/scalefactor;
output(2)=utility(rho,consu2)+beta*wint(2)/scalefactor;
output(3)=utility(rho,Y(j)-consu2-consu1)+beta*wint(3)/scalefactor;
f(1)=gammarfun(1) - consu2/consu1;
f(2)=gammarfun(2) - (Y(j)-consu2-consu1)/consu1;
f(3)= waut(j,1)-output(1);
f(4)= waut(j,2)-output(2);
end

end